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SUMMARY 


Implicit approximate factorization techniques (AF) are investigated for 
the solution of matrix equations resulting from finite-difference approxima- 
tions to the full potential equation in conservation form. For transonic 
flows, an artificial viscosity, required to maintain stability in supersonic 
regions, is introduced by an upwind bias of the density. Two implicit AF pro- 
cedures are presented and their convergence performance is compared with that 
of the standard transonic solution procedure, successive line overrelaxation 
(SLOR). Subcritical and supercritical test cases are considered. The results 
indicate that the AF schemes are substantially faster than SLOR. 


I. INTRODUCTION 


There are basically three formulations for inviscld transonic flows. 

These are, in order of increasing complexity (and in order of decreasing 
approximation): (1) transonic small disturbance potential equation (TSD) , 

(2) full potential equation (FP), and (3) Euler equations (exact inviscid for- 
mulation). TSD is valid for thin wings at free-stream Mach numbers near unity 
and is an isentropic and irrotational formulation. It offers the advantage of 
simplicity, especially in the treatment of wing surface boundary conditions. 
The FP formulation can be considered exact under the assumptions of irrota- 
tional and isentropic flow. These assumptions, which are less restrictive 
than those for TSD, are valid for a wide range of practical transonic flows. 
Potential formulations can be written in terms of a single second-order par- 
tial differential equation (PDE) , whereas the Euler formulation consists of a 
set of first order PDE's (usually four in two-dimensional cases). Hence, for 
implicit AF schemes the potential formulations require only scalar matrix 
operations, while the Euler formulation requires time— conssuming block matrix 
operations. The FP formulation is the most efficient of the three formula- 
tions in terms of accuracy-to-cost ratio for the majority of inviscid tran- 
sonic flow applications. 

The object of this investigation has been to determine the feasibility of 
using implicit approximate factorization algorithms (AF) to solve the full 
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potential equation in conservation form for steady transonic flow fields. The 
two AF schemes tested are logical extensions of the schemes previously devel- 
oped for the TSD equation (refs. 1,2). These previous studies found the AF 
approach to be substantially faster than the standard transonic flow field 
solution procedure, successive line overrelaxation (SLOR). 

Section II begins with a discussion of the full potential equation in 
conservation form. Spatial difference approximations arc then introduced, 
including the addition of artificial viscosity required to maintain stability 
in supeLSonic zones. The difference approximations are equivalent to those 
introduced by Jameson (ref. 3). Here, however, the artificial viscosity is 
not introduced explicitly, as in the Jameson approach, but rather by retarding 
the density. 

In Section III, two implicit AF Iteration schemes are presented. AF and 
SLOR convergence histories are then compared in Section IV for both subcriti- 
cf 1 and supercritical cases. Results indicate a substantial increase in com- 
putational efficiency for the AF schemes over an SLOR scheme. 

The authors express their gratitude to Joseph L. Steger for his many 
helpful suggestions during the course of this study. 


II. SPATIAL DIFFERENCING OF THE FULL POTENTIAL EQUATION IN CONSERVATION FORM 


A. The Full Potential Equation 


The full potential equation written in conservation-law form is given by 


(P*x>x + <?Vy = ° 


( 1 ) 


where 


P 



Y - 1 

Y + 1 




1/Y-l 


( 2 ) 


In equations (1) and (2) the density (p) and velocity components ($ x and $y) 
are nondimensionalized by the stagnation density (p s ) and the critical sound 
speed (a*), respectively, x and y are Cartesian coordinates, and y is the 
ratio of specific heats. In addition, other useful relations nondimensional- 
ized by p s and a^ are given by 



mi 
1 2 Y - 1 

(3) 

Y + 1 y 
P 2Y p 

(4) 

2 = IE 
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where q is the flow speed + 4>y 2 ), a is the speed of sound and p 

is the pressure. 

Equations (1) and (2) express mass conservation for flows that are isen- 
tropic and irrotational. The corresponding shock-jump conditions are valid 
approximations to the Rankine-Hugoniot relations for many transonic flow 
applications. A comparison of the isentropic and Rankine-Hugoniot shock polars 
is given in reference 4. 

It is essential that the finite difference approximation to equation (1) 
be cast in conservation form (ref. 5). Otherwise, the shock capturing pro- 
cedure will not necessarily conserve mass across the shock wave (ref. 3). 
Nonconservative rather than conservative difference schemes have been used in 
many engineering applications. However, the nonconservative procedures intro- 
duce mass sources at shock waves, and the strength of these sources depends on 
the local grid spacing, a nonphysical consideration. Erroneous shock solu- 
tions therefore result. 


B. Spatial Differencing in One Dimension 
To begin with, consider the one-dimensional version of equation (1), 

( p 4’ x ) x = 0 (6) 

A second-order accurate finite-difference approximation to equation (6) is 
given by 


$ p. . , , = 0 (7) 

x' i-t-(i/2) x T i K/J 

where and 6 X arc backward and forward difference operators defined by 

^i = Ax -1 (<ji^ - 4> i _ 1 ) 

. ( 8 ) 

Vi = Ax ’ 1 (i(, i+i - V 

Equation (7) is a suitable finite-difference scheme for subsonic flow 
regions; however, for supersonic regions, a properly-chosen artificial vis- 
cosity term must be added. For example, Jameson (ref. 3) adds the following 
viscosity term: 


where U = min[0,p(l — $ x -/a 2 )]. This is analogous to the switching used in 
the Murman mixed— difference procedure (ref. 6). It can be shown (ref. 3) (by 
differentiating the one-dimensional form of equation (2)) that this is equiva- 
lent to adding 


- (vp xVx 


( 10 ) 
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where v = max[0, (1 - a z /<t> x ? -)l* The complete finite-difference approximation 
to equation (6) can thus be wriLcen 


‘"Vi'Vwi /!) 1 


X 1 


■ Vi [P i+(l/2) 


" P i-(l/2) J Vl = 


( 11 ) 


(Other difference approximations to equation (10) are possible, and some of 
these will be discussed in a subsequent section.) This scheme is centrally 
differenced and second-order accurate in subsonic regions. In supersonic 
regions, the differencing is a combination of the second-order accurate cen- 
tral differencing used in subsonic regions and the first-order accurate upwind 
differencing resulting from the addition of artificial viscosity. As the flow 
becomes increasingly supersonic, the scheme is increasingly retarded in the 
upwind direction. 


Equation (11) can be rearranged to give 


(p4> ) 

x x 


£ o 


X M i+(l/2) x^i 


= 0 


where 


P i+(l/2) . (1 " V P i+(l/2) + Vi~(l/2) 


(12a) 


(12b) 


The addition of the artificial viscosity given by equation (10) is thus equiv- 
alent to retarding the density in equation (7). Artificial viscosity is not 
added explicitly as in the Jameson (ref. 3) procedure. However, the two 
approaches produce identical results. The significance of the present 
approach will become apparent in the discussion of solution algorithms 
(Section III). 

The choice of v strongly affects the accuracy and stability of solu- 
tions to equation (12). The particular choice v = max[(l - a 2 /4> x 2 ),0] results 
in an effective artificial viscosity that corresponds exactly to the form used 
by Jameson (ref. 3), and generally satisfactory results are obtained. For 
v = 0, the difference scheme reduces to equation (7), which is unstable for 
supersonic regions. The choice v = 1 results in a first-order-accurate 
approximation that is highly dissipative, However, stable solutions can be 
obtained for both subsonic and supersonic regions. Thus, by sacrificing 
second-order accuracy in subsonic regions, a scheme can be constructed that 
need not be switched in the Murman fashion. Moreover, accuracy can be 
improved somewhat by decreasing v. However, in supersonic regions a stabil- 
ity limit is encountered for values well above zero. Computed results for 
various choices of v are presented in the following section. 


C. Computed Results for the One— Dimensional Case 

The effect of v on shock wave resolution is studied here using a one- 
dimensional test problem illustrated in figure 1. A uniform grid spacing of 
unity is used. Boundary conditions are 
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■Ho) = o 


\ 


* x (0) = u„ > 

i|f(XMAX) = $ (specified) / 


(13) 


The slope discontinuity in the ij> distribution represents a shock wave. The 
flow upstream and downstream of the shock is uniform. 

Computed results are shown in figure 2. Three different shock strengths 
are presented, corresponding to three different choices of U^. Exact (dis- 
continuous) solutions are indicated by the dashed line. The discrete, data are 
solutions to equation (12) for three different choices of v, representing 
three different strategies for choosing a value of artificial viscosity by 
retarding the density. The re.sult equivalent to the Jameson viscosity, which 
is switched on (in the Murman fashion) only for supersonic points, is shown by 
the triangles. There is a slight overshoot that increases in amplitude with 
increasing shock strength. Note that the shock is always captured in a dis- 
tance of about three mesh cell widths. 


The two other solutions were computed with v constant. The density is 
retarded by the same amount in both subsonic and supersonic regions, so that 
there is no Murrnan-type switching. The solutions represented by the squares 
were computed using a constant value of v equivalent to that used in the 
(uniform) supersonic region by the Jameson strategy. Shock waves for this 
approach are captured in about five to seven mesh cell widths. For the choice 
v - 1, shocks are smeared over a substantially greater distance, especially 
for the weaker shocks. 

The results shown in figure 2 indicate that the Jameson choice of arti- 
ficial viscosity, implemented here by retarding the density, is a suitable one 
for shocks of weak to moderate strength. Furthermore, if one is willing to 
sacrifice second-order accuracy in subsonic regions and accept a greater 
degree of shock smearing, solutions can be obtained without switching differ- 
ences; the tesc for supersonic flow can thereby be eliminated. 


D. Spatial Differences in Two Dimensions 


A second-order act .rate finite-difference approximation to equation (1) 
is given by 


[ 6 p 




x i+(i/V),j x y i,j+(l/2) y J i,j 


ti*. 


= 0 


(14) 


where ciy and £y are backward and forward y-direction difference operators, 
respectively defined similarly to the x-direction operators given by equa- 
tion (8). As in the one-dimensional case, equation (14) is a suitable finite- 
difference scheme for subsonic flow regions but not for supersonic flow 
regions. A properly chosen artificial viscosity term must bo added. In two 
dimensions, Jameson (ref. 3) adds an artificial viscosity terra of the follow- 
ing form (u i+(l/2)>j > 0, v ifj+(l/2) > 0) 
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(15a) 




‘l.j 


y i.j 




where v = max[0,(l - a 2 /q 2 )]. For the case when vi f j+(i/ 2 ) < 0 the artifi- 
cial viscosity term is slightly different and is given by 


"^x V i,j U l+(l/2),j P x 


i.j 


~ Vi.j+l V i.j+('/2) P y i 

* 


j+1 


(15b) 


Other cases arise when u i+(i/2),j < ® and are handled in a similar fashion. 
For convenience of presentation, only the first case will be considered here- 
after. The complete finite-difference approximation to equation (1) can thus 
be written 

(»+*>* + <P Vy “ fVi+(l/ 2 ),/x + 

“ ^x v i,j u i+(l/ 2 ),j p x i ^' 6 y V i, j V i ( j+(l/2) P y 1 ^ " 0 


This scheme is centrally differenced and second-order accurate in subsonic 
regions. In supersonic regions, the differencing is a combination of the 
second-order accurate central differencing used in subsonic regions and the 
first-order accurate upwind differencing resulting from the addition of arti- 
ficial viscosity. As the flow becomes Increasingly supersonic, 
increasingly retarded in the upwind direction. 

As in the one-dimensional case, the two-dimensional scheme 
rearranged to give 

(P *A + (e Vy * 'Vi+tl/iA + Vi+d/2)Vh,j 

P i+ ( 1 / 2 ) = (1 " V i,j )P i+0/2),j + V i,j P i-(l/2),j 

P j+(l/2) = (1 " V i,j )P i,j+(l/2) + V i,j P i,j-(l/2) 

The addition of the artificial viscosity given by equation (15) 
equivalent to retarding the density in equation (14). Artificial viscosity is 
not added explicitly as in the Jameson procedure. However, the two approaches 
produce identical results. As pointed out by Jameson (ref. 3), the difference 
scheme given by equation (16) provides automatic upwind differencing of the 
streamwise terms in supersonic regions. Thus the full effect of rotated dif- 
ferencing is included in the present finite-difference scheme. 

The present scheme can easily be extended to include an arbitrary coordi- 
nate system. For example, transformation of the full potential equation from 
Cartesian coordinates (x,y) given by equation (1) to a general coordinate 
system (£,n) yields the full potential equation in the following form: 



the scheme is 
can be 

= 0 (17a) 

(17b) 
(17c) 

is thus 


* 
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where J Is Che Jacobian of the transformation and U and V are the contra' 
variant velocities along the £ and n directions, respectively, which are 
defined by 


U - 

v 

V = A,,^ + A 3 + n 


(19) 


A,, A 2 , and A<j are metric coefficients which depend on the type of transfor- 
mation. No restrictions have been placed on the (£,n) coordinate system, and, 
for example, nonorthogonal coordinates can be used. This allows considerable 
flexibility in treating a wide range of geometries. 


The spatial differencing scheme for a general (£,n) coordinate system car. 
be written as follows 


(*?),. * (*?), 


V’i+(i/2 




+ r 


i+(l/2),j 


n p j+(i/2 


>a. 


= o 


j-Kl/2) 


(20) 


where Pi+(i/2) ancl Pj+(l/2) are defined as before, by equation (17). The 
general spatial difference scheme given by equation (20) contains all of the 
properties of the simpler Cartesian version (eq, (17)). Namely, this scheme 
is second-order accurate and centrally differenced in subsonic flow regions. 
The artificial viscosity is obtained solely by an upwind bias of the density 
coefficient in both the £ and n directions. This form of spatial differ- 
encing fits nicely into the framework, of many iteration procedures. Three 
procedures will be presented and tested in subsequent sections of this report. 
For the sake of simplicity, only the Cartesian form of the full potential 
equation will be investigated. 


III. IMPLICIT APPROXIMATE FACTORIZATIONS 


A. General Requirements 

The finite-difference approximation to the full potential equation intro- 
duced in the previous section (eq. (17a)) can be written 


Li|> = 0 


( 21 ) 


and is applied to every mesh cell in the flow field. The purpose of this 
section is to describe several iteration procedures for solving the resulting 
matrix equation to determine <f> (and thus p) at each mesh point. 

A general form for a two-level solution procedure is given by 

NC n + u)R n = 0 (22) 
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where C n (=i}> n - 4> n ) is the correction, R n (=L4» n ) is the residual, which is a 
measure of how well the finite-difference equation is satisfied by the nth 
level solution (d? 11 ) , and m is a relaxation parameter. The iteration scheme 
given by equation (22) can be regarded as an iteration in pseudo-time, where 
the n superscript indicates the time-step level of the solution. 

The operator N determines the type of iterative procedure and, there- 
fore, determines the rate at which the solution procedure converges. In the 
approximate factorization approach N is chosen as a product of two or more 
factors indicated by 


N = Nj • N 2 ~ L (23) 

The factors Nj and N ? are chosen so that: (1) their product is an approxi- 

mation to L, (2) only simple matrix operations are required, and (3) the 
overall scheme is stable. 

In what follows, three iteration procedures are discussed, each corre- 
sponding to a different choice for N. These methods, which Involve only 
simple bidiagonal or tridiagonal matrix operations, are: (1) successive line 

overrelaxation (SLOR), the standard transonic flow solution procedure; 

(2) alternating direction implicit (ADI), one type of implicit approximate 
factorization scheme (called AF1 in ref. 1); and (3) AJ2, another type of 
implicit approximate factorization scheme. 


B. Successive Line Overrelaxation (SLOR) 

The SLOR algorithm used in the present study can be expressed by choosing 
N as follows 


NC? . = Ax" 1 


~n 


/ P i+(l/2) ~n \ 

\ Ax p i-(l/2) 6 xj 


Vj+a/2)‘ 


k 

y i.i 


(24) 


This difference expression has been time-linearized by evaluating p at iter- 
ation level n. Because the cross terms are indirectly included in p, the 
form of N chosen here is simpler than the commonly used quasi-linear form of 
the full potential equation. This scheme is implicit in the y direction; 
that is, the complete y-direction operator is included in N. This requires 
the inversion of a tridiagonal matrix equation for each x = constant line. 

The operator is explicit in the x direction because it contains only the 
lower diagonal part of the x-direction operator. This means that each grid 
point is influenced by only a single grid point to the right in the x direc- 
tion during one Iteration, which contributes to a relatively slow evolution of 
the solution. Construction of N in a fully implicit manner, however, means 
that each grid point is influenced by every other grid point during each 
iteration. As a result, much faster convergence can be obtained. The next 
two algorithms presented are both fully implicit and obtain this form by con- 
structing N as approximate factorizations of L. 
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C. Approximate Factorization, Scheme 1 (AF1) 

The first fully implicit algorithm is similar to the AF1 scheme presented 
in reference 1, which was used to solve the TSD equation. This scheme, which 
is a reformulation of the Peaceman-Rachf ord alternating direction implicit 
technique (ADI), can be expressed by choosing N as follows: 


uNC? . 
i.J 



- t 


,.n 

x p i+(l/2) 




j+O/'O 



n 

i.j 


(25) 


where a is a free parameter to be defined subsequently. Note that both the 
x and y directions are treated implicitly and that N has been written as 
the product of two factors which when multiplied out yield 


i,j i ,3 x l+(l/2) x y ji(l/2.) y x,3 


+ aLC 


n 




(26) 


This expression includes the time-linearized L operator plus two error terms. 
The first error term is a t — type term and therefore, provides a stabilizing 
effect to the iteration process in subsonic flow regions but a destabilizing 
effect in supersonic flow regions (ref. 7). The effect of this term is appar- 
ent in the computed results, Section IV. 


The scheme can be restated in practical terms using two steps as follows 


step 1 : 


[« 



1 + 0 / 2 ) 



n 

i» j 


oui)L<}> 


n 

i.j 


(27a) 


step 2 ; 


[« - 6 P 1 ?./, ]C? . = f? . 
y 0+0/2) y J X,J 


(27b) 


where to is a relaxation parameter and f^ j is an intermediate result 
stored at each mesh point in the finite difference mesh. In step 1 the f 
array is obtained by solving a tridiagonal matrix equation for each y = con- 
stant line. The correction array (cj s j) is then obtained in the second step 
from the f array by solving a tridiagonal matrix equation for each x = con- 
stant line. Thus by writing N as the product of two factors it is possible 
to obtain a fully implicit technique involving only simple tridiagonal matrix 
operations. 


D. Approximate Factorization, Scheme 2 (AF2) 


The second fully implicit algorithm presented here Is similar to the AF2 
scheme presented in reference 1, which was applied to the TSD equation. The 
AF2 scheme, first investigated by Ballhaus and Steger (ref. 8) for unsteady 
flows, can be expressed by choosing N as follows: 


aNC? . 



- Vi+o/o^X 0 


5 i+(l/2) 




(28) 
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where again N has been written as the product of two factors, which, when 
multiplied out, >ield the time-line;.; rized L operator plus two error terms. 
The first error term is an upwind <}> xt -type term and therefore, provides time- 
dependent dissipation for the convergence process, which is especially conven- 
ient for supersonic flow. When constructing N in this form, care must be 
taken to ensure that this term has the proper sign and is differenced in the 
upwind direction. Implementatic n of the AF2 scheme is achieved by writing it 
in a two step form given by 


step 1 : 


step 2 : 


t Jf" . 
y/ i.j 

= . 

3 

(29a) 

*y? < 

> i.j 

= £ n 
±,i 

(29b) 


where to is a relaxation parameter and ^i,j is an intermediate result 
stored at each mesh point in the finite difference mesh. In step 1, the f 
array is obtained bv solving a tridiagonal matrix equation for each x - con- 
stant line. The correction array (C^j) is then obtained in the second step 
from the f array by solving a simple bidiagonal matrix equation for each 
y = constant line. Note that with AF2 the x-direction difference approxima- 
tion is split between the two steps. This generates the desired $ X £ term as 
mentioned above and also places a sweep direction restriction on both steps, 
namely, downwind for the first step and upwind for the second step. 


Several variations of the AF2 scheme are possible. Some of these fnclude: 

(1) splitting the y-direction term instead of the x-direction term, 

(2) changing the order of the steps, and (3) moving the x-direction density 
coefficient (Pi-f(i/ 2 )) from the second factor to the first factor. These 
three variations ha^e been tried and were found to be stable. 


Normally, flow-field, type-dependent differencing is used to achieve 
stability in transonic flow calculations. Incorporating these different oper- 
ators into iteration procedures, such as the AF schemes presented here, would 
be cumbersome if not impossible. Using the upwind bias of the density coeffi- 
cient, which is always evaluated at the nth iteration level, allows the 
simple two- and three-banded matrix form of the AF schemes to be retained over 
the entire flow field, even in regions of supersonic flow. In fact, use of 
upwinded density coefficients in any general iteration scheme (i.e., for any 
arbitrary N operator) could be used to remove the difficulties introduced by 
type-dependent differencing. The resulting general scheme would retain the 
same basic differencing (at the n+1 iteration level) throughout the entire 
flow field, relying on the upwind bias of the density (at the nth iteration 
level) to provide the artificial viscosity in supersonic flow regions. This 
represents a significant simplification in the handling of supersonic flow 
regions for transonic flow calculations. 


E. a Selection for the Approximate Factorization Schemes 

In both the AF1 and AF2 schemes the quantity a is an as yet undefined 
free parameter. If a were chosen to be At -1 then both schemes could be 
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considered to be iterations in pseudo time. This provides one strategy for 
obtaining fast convergence, namely, advance time as fast as possible with 
large time steps (i.e., small ct's). As pointed out in reference 1, this is 
effective for attacking the low-frequency errors but not the high-frequency 
errors. The best overall approach is to use an a sequence containing several 
values of a. The small values are particularly effective for reducing the 
low-f requeue'/ errors, and the large values are particularly effective for 
reducing the high-frequency errors. Endpoints for a suitable u sequence can 
be approximated analytically; for the present case these approximations are 
given by 


AF1 : 

a H = 4p/Ay 2 

a L - P 

(30) 

AF2 : 

a H = 1/Ay 

« L = 1 

(31) 


where p is a representative value of density (e.g., p^) , and Ay is the 
minimum y-direction spacing in the finite-difference mesh. Refinement of 
these estimates by numerical experiment increases the efficiency of the itera- 
tion procedure. The a sequence used in the present study is given by 

(k-T )/ (M-l) 

| k = 1, 2, . . . , M (32) 

if*, 

where M is the number of elements in the sequence. This (geometric) 
sequence has been found to be effective (ref. 1), although other sequences 
could perhaps provide equivalent or improved convergence performance. 



IV. TWO-DIMENSIONAL RESULTS 


The three schemes presented in the previous section (SLOR, ATI, and AF2) 
are evaluated in this section. A two-dimensional, 10% thick circular-arc 
airfoil with small-disturbance boundary conditions is used as a test case. 

Both subcritical (Case A) and supercritical (Case B) Mec n numbers are consid- 
ered (see table 1). The finite-difference grid is Cartesian with variable 
spacing in both the x and y directions, as shown in figure 3. Use of the 
Cartesian grid with small-disturbance boundary conditions was motivated from 
the standpoint of simplicity and does not reflect the desired ultimate use of 
the full potential schemes under investigation. 

The subcritical and supercritical pressure coefficient distributions 
(Cases A and B, respectively) are presented in figure 4. These results were 
computed with a 90 * 21 mesh containing 47 points on the airfoil surface. The 
boundaries were located 5 chord-lengths away from the airfoil in the x 
direction and 6 chord-lengths away in the y direction. In the AE cases, the 
a sequence contained eight elements (M=8) . 
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A. Convergence Performance for the Subcritical Case 


Convergence characteristics for the subcritical case in terms of maximum 
residual histories are displayed in figure 5. The AF results were established 
by plotting the maximum residual every eighth iteration, which always corre- 
sponded to a = (xl in the eight element sequence. The SLOR results were 
established by plotting the residual every 20 iterations. All of the conver- 
gence parameters (mgL0R> a H» anc * a L^ h ave been selected by a trial and error 
optimization process. (Results for two choices of relaxation parameter 
(^SLOR ~ 1-95 and 1.975) are presented for SLOR. The larger value produced 
faster convergence for a maximum residual drop greater than four orders of 
magnitude, while the smaller value produced faster convergence for smaller 
drops in the residual.) The relaxation factor used in both the AF1 and AF2 
schemes (w in eqs. (27) and (29)) was held fixed at 2 for all test cases. 

Based on a six order of magnitude drop in maximum residual, the AF1 scheme is 
about twice as fast as the AF2 scheme, and about 16 times faster than the SLOR 
scheme. These speed ratios are in terms cf iteration count. The AF1 and AF2 
schemes take about 50% and 30% more CPU time per iteration than SLOR. This 
should be taken into account when considering the speed ratios based on the 
total amount of computational work. 

For several reasons, exact determination of speed ratios for these 
schemes is difficult to assess. Firstly, use of grid sequences usually pro- 
vides as much as a factor of 2 speed increase for SLOR schemes, but only a 
small speed increase for the AF schemes (ref. 1). Grid sequences were not 
used for any results reported here. Secondly, use of nonoptimal convergence 
parameters does slow convergence by as much as a factor of 2 or more for both 
the SLOR and AF schemes. Because the AF schemes have two parameters to opti- 
mize (ay and a^) , as opposed to only one for SLOR (mgLOR^ • optimization is more 
difficult for the AF schemes. Finally, and most importantly, the AF schemes 
reduce the errors associated with all frequencies equally well, approximately, 
while the SLOR scheme is efficient for only the high frequencies. Since the 
residual is heavily biased toward the high frequency end of the error spectrum, 
using a specified drop in the maximum residual to define convergence heavily 
favors the SLOR scheme (ref, 2). More discussion on this last point is pro- 
vided in this section, part C. 


B. Convergence Performance for the Supercritical Case 

Convergence characteristics for the supercritical case (Case B) in the 
form of maximum residual histories are displayed in figure 6. Again the 
w SL0R» “H« anc * a L va lues have been obtained by trial and error optimization. 
Based on a six order of magnitude drop in the maximum residual and in terms of 
iteration count, AF2 is slightly more than twice as fast as AF1, and about 
11 times faster than SLOR. 

The number of supersonic points (NSP) plotted versus iteration number for 
Case B is shown in figure 7. The final MSP is 187. The AF2, AF1, and SLOR 
results reach this level in 29, 103, and 320 iterations, respectively. The 
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small number of Iterations required by AF2 to establish the supersonic zone is 
another indication of how rapidly the solution evolves with this approach. 

The AF2 scheme was relatively consistent In convergence speed for both 
the subsonic and supersonic cases (79 and 118 iterations, respectively). The 
AF1 scheme, on the other hand, displayed remarkable speed for the subsonic 
case (39 iterations) but was a disappointment for the supersonic case 
(254 iterations). Perhaps this is because the $ xt error term produced by 
the AF2 factorization is more suitable for supersonic regions than the <f> t 
error terra resulting from the AF1 factorization. 

It should be pointed out that use of the standard definition for v in 
two dimensions (see eq. (10)) produced pre-shock, overshoots which often 
resulted in numerical instability. This instability was experienced for both 
the standard rotated (see eq. (17)) and nonrotated difference schemes. (The 
nonrotated difference scheme simply has no y-directlon artificial viscosity 
term. This causes the term (p 4>y ) y to remain centrally differenced in super- 
sonic regions.) Therefore, an alternative definition for v y.'.s introduced, 



(33) 


where p* is the sonic value of density and a e has been chosen by numerical 
experiment to be six. Use of equation (33) instead of the standard definition 
for v increases the amount of upwinding or, equivalently, the amount of 
artificial viscosity in supersonic flow regions. Shock wave overshoots were 
thereby prevented for the supercritical case presented here. 


Several variations of the artificial viscosity term have been investi- 
gated. In all cases, only the x-direction artificial viscosity term has been 
included (i.e., nonrotated differencing). Two of these variations are given 
by 



j 


u. .p 
1,3 


(34) 


and 


^x U i+(l/2) ,j U i+(l/2) , j P x. . < - 35 

1 > J 

The artificial viscosity term introduced previously (eq. (15)) corresponds 
exactly to the term reported by Jameson (ref. 3). This form is incorrectly 
reported in that Jameson actually uses the artificial viscosity given by 
equation (34) (private communication from A. Jameson, Courant Institute of 
Mathematical Sciences, N.Y.). The artificial viscosity term given by equa- 
tion (35) is still another successful version. All three variations of arti- 
ficial viscosity (with suitably tailored forms for v) have been tested for 
Case B and give essentially the same results. 
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C. Residual vs Error 


Relative levels of convergence for AF2 and SLOR for given reductions in 
maximum residual are compared in figure 8. The solid line represents the 
final solution, which has been converged until the maximum residual dropped 
six orders of magnitude. The other results represent intermediate AF2 and 
SLOR solutions in which the maximum residual has been reduced by one, two, and 
three orders of magnitude. It is immediately obvious that reducing the maxi- 
mum residual by equal amounts for AF2 and SLOR does not produce intermediate 
results with the same level of error. In fact, the error reduction for these 
two schemes can be substantially different for the same degree of residual 
reduction. For instance, the AF2 solution after a two order of magnitude 
reduction in maximum residual (fig. 8(b)) is closer to the final solution than 
the SLOR solution is after a residual droD of three orders of magnitude 
(fig. 8(c)). 


This behavior can also be observed by comparing the maximum residual his- 
tory curves of figure 6 with the RMS error history curves given in figure 9. 
The RMS errcr at iteration n (Ef^g) was computed from the surface pressure 
coefficient distribution by the following formula 



(36) 


were C is the surface pressure coefficient at the ith grid point and the 
H i 

nth iteration, C p is the surface pressure coefficient at the ith grid 


point taken from the converged solution, and i^g and iyg are the indices 
indicating the leading and trailing edge grid points. The SLOR residual drops 
very rapidly initially and then levels off. The SLOR RMS error drops grad- 
ually. Therefore, at the ’’knee" in the residual history curve, even though 
the residual has dropped by about three orders of magnitude, the actual RMS 
error has dropped by only one order of magnitude. In contrast, both maximum 
/residual and RMS error results for the AF schemes are nearly straight lines 
with about the same slope. 


This behavior is the result of two factors (ref. 2): (1) the AF2 scheme 

treats all error components equally well (approximately) , whereas the SLOR 
scheme performs efficiently on only the high-frequency error components; and 
(2) it can be shown that the residual is a weighted sum of errors, over the 
entire error frequency spectrum, weighted by the eigenvalue of the finite 
difference scheme. The eigenvalue for the high-frequency error components is 
0(Ax 2 ), while the eigenvalue for the low-frequency error components is 0(1). 
Hence, the residual is heavily influenced by the high-frequency errors. During 
the initial phase of a calculation, equal residual drops for AF2 and SLOR 
indicate the same drop in high frequency error, but the reduction rates for 
the low-frequency errors are significantly different. For this reason the AF 
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approach reduces the total error much faster, and consequently, the solution 
evolves much more rapidly. Therefore, the maximum residual operator should 
not be used as the basis for comparison between AF and SLOR schemes. RMS 
errors are much better suited for this purpose. In practice, use of the maxi- 
mum residual operator for either AF or SLOR is the most convenient method for 
monitoring convergence (since the error is unknown). However, the convergence 
criterion must be adjusted in accordance with the solution procedure in use. 


D. Use of Nonoptimal Convergence Parameters 

Only "optimal" convergence (l.e., cases with convergence parameters 
adjusted for optimal convergence) have been considered thus far. Because in 
practice the optimal values for ay, a^, and are not known a priori, it 

is of interest to know the effect of nonoptimal parameters on the convergence 
speed. The RMS error histories for AF1, AF2, and SLOR for nonoptimal param- 
eters are shown in figure 10. The nonoptimal SLOR relaxation factor was 
chosen to be 1.8, and the and parameters for AF1 and AF2 were chosen 
from equations (30) and (31). For a four order of magnitude reduction in RMS 
error (which is well beyond plot table accuracy) nonoptimal AF2 is approxi- 
mately 2.5 times slower than the optimal AF2 results. The AF1 and SLOR results 
are similarly affected, being about 1.5 and 3 times slower, respectively. For 
the nonoptimal case, in terms of iteration count, AF2 is about 1.5 times 
faster than AF1 and about 12 times faster than SLOR. In terms of actual com- 
puter time, AF2 is about 1.7 times faster than AF1 and nearly an order of 
magnitude faster than SLOR. 

Another indication of how fast the AF2 scheme can establish the global 
solution even for nonoptimal acceleration parameters is shown in figure 11. 

The solid line in each case is the final solution, which has been converged 
until the maximum residual dropped six orders of magnitude. The other results 
presented are AF2 and SLOR intermediate solutions in which the maximum resid- 
ual has been reduced by one, two, and three orders of magnitude. As before 
(fig. 8), it is again obvious that reducing the maximum residual by equal 
amounts for AF2 and SLOR solutions does not produce equivalent error reduction. 
In the present nonoptimal case, the difference between the AF2 and SLOR solu- 
tions for a given drop in residual is even larger than in the previous optimum 
case. This is because the smaller value of mg^QR used in the nonoptimum 
case puts even more emphasis on the high-frequency errors, thus causing the 
SLOR residual to drop faster, while the RMS error is dropping slower. 

Obtaining nearly optimal a-sequence endpoints (a^a^) for the AF2 scheme 
is not difficult. The most important influence on the a values is the 
finite— difference mesh. If the mesh remains fixed from one case to another, 
then the same a's will still be approximately optimum. If the mesh does 
change, estimates for the new optimal a endpoints (a£ eW , afp”) can be 
obtained (from eq. (31)) as follows: 
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( 37 ) 


new _ old 


new old . old,, new 
a = i Ay . /Ay . 

H H ^min min 

where and ‘h^ are tlie °ld hand-optimized (trial and error) endpoints, 

and Ay^ in and AyJ^ are fh e minimum Ay spacings in the old and new meshes, 
respectively. This philosophy has been tried in a limited number of cases, 
and seems to be acceptable. However, because of the limitations of the pres- 
ent formulation (nonlifting and small disturbance boundary conditions) more 
investigation is needed. 


V. CONCLUSIONS 


New fully implicit algorithms for solving the full potential aquation in 
conservation form have been developed. These new schemes are of the 
approximate-factorization variety. Computed results indicate a substantial 
increase in computational efficiency over an SLOR algorithm. 

The spatial difference scheme used to approximate the full potential 
equation is equivalent to the one developed by Jameson. However, in the pres- 
ent approach, the artificial viscosity required to maintain stability in 
supersonic regions is not added explicitly. It is introduced by spatially 
retarding the density in the finite-difference equation. This strategy 
greatly simplifies the solution procedure, so that only bidlagonal or tridiag- 
onal matrix operations are required. 

Results indicate that the standard measure of convergence (i.e., a speci- 
fied reduction in the maximum residual) is not a good means for comparing AF 
and SLOR convergence rates. Using the maximum residual as a criterion to 
determine relative levels of convergence for either technique individually, is 
not questioned. However, AF schemes definitely reach a greater degree of con- 
vergence in terms of the error at higher residuals than do SLOR techniques. 

The new AF schemes although applied in this report only to circular-arc 
airfoils with small-disturbance boundary conditions and Cartesian finite- 
difference grids, are extendible to arbitrary geometries in both two- and 
three-dimensional flows. The principal difficulty is the handling of the 
cross derivative terms arising from general geometry transformations. Solu- 
tions for cases with airfoil adapted grids have been obtained, and the effi- 
ciency of the AF procedure does carry over to these more practical flows. 
Details are presented in a report in preparation by T. L. Holst. 
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TABLE 1.- TABULATION OF SOLUTION PARAMETERS 



Case a 

Scheme 

“SLOR 

“L 

a H 



(SLOR 

1.975 


— — — 

Optimum 

A 

{afi 

— 

0.04 

100,000 



1 AF2 

— 

.4 

100 



(SLOR 

1.95 

— — 

— 

Optimum 

B 

{afi 

— 

1.3 

4,000 



|aF2 

— 

.6 

60 



( SLOR 

1.8 


— 

Nonoptimum 

B 

{afi 

— 

.8 

128,000 



{ AF2 

— 

1.0 

200 


a Case A: M^O.? (subcritical) ; Case B: M m = 0.84 

(supercritical) . 






Figure 2.- Shock profiles, <j> x vs x. 
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Figure 3.- Finite-difference mesh (two dimensions). 
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Concluded. 
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Figure 5.- Maximum residual conv>r? 1 enre history comparison (Case A). 



Figure 6.- Maximum residual convergence history comparison (Case B) 
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Figure 8.- Intermediate solution comparisons after specified reductions in the 
maximum residual (Case B, optimum convergence). 
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Continued. 



Concluded 



30 


Figure 9.- RMS error convergence history comparison (Case B). 
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Figure 11.- Intermediate solution comparisons after specified reductions in 
the maximum residual (Case B, nonoptimal convergence). 




Figure 11.- Continued. 
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